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ABSTRACT 

We propose a sampling scheme on the sphere and develop a cor¬ 
responding spherical harmonic transform (SHT) for the accurate re¬ 
construction of the diffusion signal in diffusion magnetic resonance 
imaging (dMRI). By exploiting the antipodal symmetry, we design 
a sampling scheme that requires the optimal number of samples on 
the sphere, equal to the degrees of freedom required to represent 
the antipodally symmetric band-limited diffusion signal in the spec¬ 
tral (spherical harmonic) domain. Compared with existing sampling 
schemes on the sphere that allow for the accurate reconstruction of 
the diffusion signal, the proposed sampling scheme reduces the num¬ 
ber of samples required by a factor of two or more. We analyse the 
numerical accuracy of the proposed SHT and show through experi¬ 
ments that the proposed sampling allows for the accurate and rota- 
tionally invariant computation of the SHT to near machine precision 
accuracy. 

Index Terms — diffusion magnetic resonance imaging; sam¬ 
pling; spherical harmonic transform; antipodal signal; unit sphere. 

1. INTRODUCTION 

Diffusion magnetic resonance imaging (dMRI) allows for the struc¬ 
ture and connectivity of white matter in the brain to be determined 
non-invasively for the detection of neurological diseases and in pre¬ 
operative planning Q][2j. White matter has a fibrous structure; the 
myelin sheath around fibre axons hinders the movement of water 
molecules, constraining diffusion. The diffusion of water molecules 
in white matter can therefore be used to determine the location and 
orientation of fibres. 

Acquisition and reconstruction of the diffusion signal from q- 
space measurements, where q is the diffusion wave vector, is a cen¬ 
tral problem in dMRI. It is common for samples to be obtained on 
a single sphere or multiple spheres in q-space (3}{5|. The number 
of measurements in q-space (images) that can be acquired in dMRI 
is severely limited due to the need for the scan time to be practical 
in a clinical setting. Hence, a q-space sampling scheme that allows 
for the accurate reconstruction of the diffusion signal using the min¬ 
imum number of samples is desirable j6|7j. 

The reconstruction of the diffusion signal can be carried out by 
expanding the signal in terms of spherical harmonics |8}|10), which 
serve as natural basis functions for signals defined on the sphere ED- 
By choosing a sufficiently large band-limit, the diffusion signal can 
be represented in terms of a finite number of spherical harmonic co¬ 
efficients, which are obtained using the spherical harmonic trans¬ 
form (SHT) calculated numerically from a finite number of measure¬ 
ments of the diffusion signal obtained at a constant q-space radius 
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(on a sphere) 0U3- In order for the accurate reconstruction of the 
diffusion signal, a sampling scheme on the sphere for obtaining q- 
space measurements needs to allow for the accurate computation of 
the SHT. 

1.1. Relation to Prior Work 

Many sampling schemes on the sphere in dMRI use an antipodally 
symmetric sampling grid to reduce the number of measurements; 
the antipodal symmetry of the diffusion signal enables the signal to 
be evaluated on the sphere at locations antipodal to where measure¬ 
ments have been obtained sun). Most of these schemes are de¬ 
signed to obtain uniform (near equal area) sampling on the sphere to 
ensure that the reconstruction is rotationally invariant, which means 
that the accuracy does not vary significantly if the sampling grid (or 
diffusion signal) is rotated (3||l7). Other sampling scheme designs 
that focus on the accurate computation of the SHT rather than geo¬ 
metric properties of the diffusion signal have also been used (5). 

In order to acquire samples of the diffusion signal band-limited 
at L (formally defined in Sectio n|2. 1 f over multiple spheres in q- 
space, an equiangular scheme ||18| that requires asymptotically 
2 L 2 samples to accurately compute the SHT has been used |5|. 
In comparison to the geometric schemes, this scheme 118] has an 
iso-latitude arrangement of samples (where samples along longitude 
are taken over iso-latitude annuli) that enables the fast computation 
of the SHT. However, the samples in the scheme are non-uniform 
due to dense sampling near the poles, therefore it is considered that 
the accuracy of the SHT changes if the signal is rotated ]3| |17|. The 
spherical design with uniform density sampling method in ED has 
a uniform and antipodally symmetric arrangement of samples, and 
allows for the accurate computation of the SHT using least squares. 
The number of samples required by this scheme does not have a 
general formula depending on the band-limit but as reported in [19] 
is greater than L~. 

For the accurate computation of the SHT, the minimum num¬ 
ber of samples, referred as the optimal spatial dimensionality, re¬ 
quired by any sampling scheme is equal to the number of spheri¬ 
cal harmonic coefficients required to represent a band-limited sig¬ 
nal in spherical harmonics basis | |18[|20| . An iso-latitude scheme 
that attains optimal spatial dimensionality ( L 2 samples) of an ar¬ 
bitrary band-limited signal on the sphere has been developed (20| , 
which allows the accurate and fast computation of the SHT. How¬ 
ever, the samples in the sampling scheme are neither antipodally 
symmetric nor uniform (equal areal) by design. Due to antipodal 
symmetry of the diffusion signal, the band-limited diffusion sig¬ 
nal has No = L(L + l)/2 degrees of freedom (see Section[3|, 
which is referred to as the optimal spatial dimensionality of the dif¬ 
fusion (antipodal) signal. We have established that the existing sam¬ 
pling schemes require at least L 2 « 2No samples on the sphere for 








the accurate computation of the SHT (or the accurate reconstruction) 
of the diffusion signal. 

1.2. Contributions 

In this work, we address whether it is possible to develop a sampling 
scheme on the sphere (at a constant q-space radius) for measuring 
the diffusion signal that 1) attains the optimal spatial dimensional¬ 
ity, i.e., requires No samples, 2) allows for the accurate computa¬ 
tion of the SHT up to and beyond practically relevant band-limits, 
3) permits rotationally invariant reconstruction, and 4) enables fast 
computation of the SHT with lowest known complexity. In address¬ 
ing these questions, we develop a sampling scheme on the sphere 
and corresponding SHT that, in contrast to previous schemes, attains 
the optimal spatial dimensionality for the accurate reconstruction of 
the band-limited diffusion signal at a constant q-space radius (an 
antipodally symmetric signal). The proposed scheme is antipodal, 
but unlike existing schemes that only focus on the geometric prop¬ 
erties of the diffusion signal, the SHT associated with the proposed 
scheme is computationally efficient and accurate, with reconstruc¬ 
tion error near machine precision (1CP 14 ) for the band-limits of in¬ 
terest in dMRI. Though the scheme is non-uniform by design, we 
show that the numerical accuracy of the scheme is the same when 
the diffusion signal on the sphere is rotated, demonstrating that the 
scheme allows for the rotationally invariant reconstruction of the dif¬ 
fusion signal. 

This paper is organised as follows. The mathematical prelimi¬ 
naries necessary to understand this work are contained in Section]?] 
In Section[3] we state the problem under consideration. The sam¬ 
pling scheme and the SHT for measurement and reconstruction of 
the diffusion signal are presented in Section]?] This scheme is then 
evaluated in terms of numerical accuracy and rotational invariance 
in Section]?] Finally, concluding remarks are made in Section]?] 

2. MATHEMATICAL PRELIMINARIES 

2.1. Signals on the Sphere and Spherical Harmonics 

Square integrable complex functions on the unit sphere § 2 have the 
form f(9, (j>), where the angles co-latitude 9 £ [0,7r] and longitude 
4> € [0, 27r) parameterise a point (sin 9 cos 4>, sin#sincos#)' £ 
R' 1 on § 2 . The inner product of two functions /, g defined on § 2 is 
given by m 

(f,g)=[ smOdOdc/), (1) 

J s 2 

where (•) denotes the complex conjugate and sin 9 d9 d(j> is the dif¬ 
ferential area element on the sphere. With the inner product in 0- 
the set of square integrable complex valued functions on § 2 forms a 
Hilbert space, denoted by L 2 (§ 2 ). The inner product in 0 induces 
a norm ||/|| = (/, f) 1 ^ 2 . We refer to functions with finite induced 
norm as signals on the sphere. 

The spherical harmonic functions (spherical harmonics for 
short) Y e m {9, <j>), defined for integer degree l > 0 and integer order 
|m| < £, form a complete orthonormal set of basis functions (HHD 
The signal is said to be band-limited at degree L if the signal can 
be completely represented in spherical harmonics basis functions 
Y e m (9, cj>) with f < L and \m\ < t. 

2.2. Diffusion Signal on the Sphere 

The diffusion signal at a fixed q-space radius can be represented 
as a band-limited signal on the sphere |12||l7| . Furthermore, the 
diffusion signal, denote by d(9, cj>), is antipodally symmetric, with 


d(9, (j>) = d{it — 9, <j> + 7r). Since Y™(9, (j>) = Y/ n (n — 9, n + (j>) 
for even i and Y[ rl (9, cj >) = —Y™(n — 9, n + cj>) for odd l, the 
expansion of d(9 , (/>) in spherical harmonic basis only includes even 
degree spherical harmonics, that is, 

L-l t 

d(9,<!>)=J2 £ ( d )T Y e m (e,<P), L odd, (2) 

£=0 m=—£ 
i even 

where ( d)™ denotes the spherical harmonic coefficient of degree i 
and order m, and is calculated using the SHT, given by 

(d)T = [ d(9, if))Y/ n (9, cf) sin 0 d9 d<j>. (3) 

is 2 

The spherical harmonic coefficients (d)™ form the spectral domain 
representation of d(9, </>). The reconstruction of the signal d(9 , (f >) 
from its spherical harmonic coefficients as given in 0 is referred as 
the inverse SHT. The band-limit required to accurately represent the 
diffusion signal depends on the q-space radius J5), and band-limits 
up to L = 11 are typically used m- 

3. PROBLEM FORMULATION 

The diffusion signal d(9,(f>) has only even order spherical harmonic 
coefficients due to its antipodal symmetry. The number of spheri¬ 
cal harmonics required to represent d(9,<j>) in 0 is No = L(L + 
l)/2|7p2), which also represents the optimal dimensionality attain¬ 
able by any sampling scheme that allows the accurate computation 
of the SHT of the diffusion signal. We customise the recently de¬ 
veloped sampling scheme p0), which requires the optimal number 
of samples for the accurate computation of the SHT of an arbitrary 
band-limited signal (without antipodal symmetry), for the acquisi¬ 
tion and reconstruction of the diffusion signal in dMRI. 

3.1. Optimal Dimensionality Sampling Scheme 

The optimal dimensionality sampling scheme in |20| has an iso¬ 
latitude sampling grid, with samples taken over L iso-latitude rings. 
The L locations along 9 where the iso-latitude rings are placed are 
stored in the vector 6 defined as 

[00, 6 U 9 L - 1 ] t . (4) 

For the accurate computation of the SHT, there needs to be at least 
2n + 1 samples along longitude in the ring placed at 9 n £ 6. The 
iso-latitude rings are placed along 9 such that the SHT can be com¬ 
puted accurately (see [20| | for further details). For the accurate com¬ 
putation of the SHT, the number of samples required by the optimal 
spatial dimensionality scheme is L 2 , which is equal to the number 
of degrees of freedom required to represent an arbitrary band-limited 
signal defined on the sphere. 

3.2. Research Questions 

In this work, we address the following questions: 

• How can we customise the structure of the optimal dimen¬ 
sionality sampling scheme to make it suitable for the acquisi¬ 
tion of measurements of the antipodal diffusion signal? 

• Can we exploit the antipodal symmetry property of the dif¬ 
fusion signal to develop a scheme that requires only No = 
L(L + l)/2 rather than L 2 samples while still allowing for 
the accurate computation of the SHT ? 










4. PROPOSED SAMPLING SCHEME 


4.1. Proposed Sampling Scheme — Structure and Design 

The diffusion signal on the sphere is antipodally symmetric, with 
d(9, <j >) = d( 7r — 9, <j> + 7r). Hence, by measuring the signal at a 
location (9a, 4 >a), we also know the value of the signal at a second 
location ( 9b,4>b) = (it — 9 a, 4>a + tr). We use this property to 
design the sampling scheme in order to reduce the number of mea¬ 
surements of the diffusion signal that need to be taken. 

With this consideration, we propose to place the iso-latitude 
rings in pairs such that they are antipodal to one another. The vector 
G, given in 0- describing the location of the iso-latitude rings is 
6 = [0,..., 7T - 9 L - 3 , 8 L - 3 , 7T - Ol- 1 ,9l-i] T , L odd. (5) 
We shortly present the location of these samples along co-latitude. 
Since we require at least 2n + 1 samples along <j> in the iso-latitude 
ring placed at 9 n for the accurate computation of the SHT, we pro¬ 
pose equiangular sampling along longitude, with fc-th sample loca¬ 
tion, denoted by </>£, in the ring placed at 8 n is given by 

n = 0 , 2, k £ [0, 2n], 

{^W 1 ’ n = l, 3, ..., L-2, k £ [0, 2(n + 1)]. 

(6) 

Remark 1. The samples in the proposed scheme are structured in a 
way that the samples in the ring 8 n are antipodal to the samples in 
the ring 8 n -\for n = 2, 4, L — 1. Therefore, we only need to 
take the measurements over the rings 9 n for n = 0, 2, ..., L — 1. 
The values of the signal over the remaining samples can be deter¬ 
mined by using the antipodal symmetry of the diffusion signal and 
exploiting the structure of the proposed sampling scheme. As an ex¬ 
ample, Fig^l^shows the proposed sampling scheme for L = 7. 

Remark 2. Since the measurements are only required to be taken 
over (L + l)/2 rings in the proposed scheme, the total number of 
measurements required by the proposed sampling scheme is 

J2(2n+1)= {L+ 2 1)L = N 0 , (7) 

n =0 
neven 

which also represent the degrees of freedom required to represent 
band-limited antipodal signal. Thus, the proposed scheme attains 
the optimal spatial dimensionality. 

4.2. Spherical Harmonic Transform (SHT) — Formulation 

We here present a variant of the method developed in |20) , to 
compute the SHT of the diffusion signal sampled using the pro¬ 
posed sampling scheme. We define an indexed vector G r " = 
[#| m |, #| m | +1 , ..., 9l-i] t C 6 for order \m\ < L, which contains 
the last L — \m\ points in the vector 6 . Also define a vector 

g m = G m (6> m ) 4 [G m (0 M ), G m (0 M+ i), ..., Gm( 6 L -i)] T , 

with 

G m (9n)= [~ W f(9 n ,<j>)e~ irn *d<p = 27T (f)T?r( 8 n), (8) 

•'° l=\m\ 

for each order \m\ < L and each 8 n £ 6 , where P( n ( 8 rl ) = 
Y™(0n, 0) denotes scaled associated Legendre functions. 

If Gm(9 n ) for each 9„ £ 6 m is computed, the spherical har¬ 
monic coefficients of order m can be recovered from 0 by setting 
up a system of linear equations |20| , given by, 

gm = P?fm, \m\ < L, (9) 



(a) (b) 


Fig. 1: (a) North pole view and (b) South pole view of the proposed 
sampling scheme on the sphere given by 0 and 0 for measuring 
d(9, (jf) band-limited at L = 7. The locations where measurements 
need to be taken are shown in black and the locations where d(9, f) 
is evaluated using antipodal symmetry are shown in grey. 


where P^ 1 is a matrix of size (L — m) x (L — m) defined as 

/ £|h(0|h) P&\+M™\) ■■■ Pl- i(*m) \ 

pr A 2?r jfa|(*IH+i) ^IH+i(*IH+i) PL-i(.0\ m \ +1 ) 


\P(Z |(0i-i) PF-i(.0l-i)J 

and f m is a vector composed of spherical harmonic coefficients of 
order \m\ < L, given by 


fm= [(/)W, (/)W+1, - •,(/)?-!] 5 


( 10 ) 


4.3. Spherical Harmonic Transform — Computation 


The spherical harmonic coefficients of order \m\ < L contained in a 
vector f m can be computed accurately by inverting 0, provided G m 
is chosen such that P'f is well-conditioned and Gm(9 n ) for each 
9 n £ 6 m is computed accurately. Since we have considered in the 
design of the proposed scheme that the ring placed at 9 n contains at 
least 2 n+ 1 samples points along longitude, Gm(9 n ) for each 9 n £ 
6 m can be exactly computed by evaluating the integral equation in 
0 as a summation over 2 n + 1 samples along longitude |20| . 

We also require the sampling points along co-latitude to be cho¬ 
sen such that the matrix PJ 1 is well-conditioned for each \m\ < L. 
Since the locations of the iso-latitude rings along co-latitude in the 
vector 9 are required to be in pairs, as given in 0. we are only re¬ 
quired to obtain measurements of the diffusion signal at (L + l)/2 
positions ( 8 q, 82 , ■ ■ ■, 9l- 1 ) along co-latitude. We define a set of 
equiangular (L + l)/2 samples along co-latitude given by 


0 = 


M2t + l) j 


t = 0,1,.. 


L — 1 
2 


( 11 ) 


and propose the following method to determine the optimal ordering 
of samples in the vector 6 : 


• Choose 9c-i = farthest from the poles, which 

is a natural choice for the ring of 2L — 1 (the largest number 
of) samples along cf>. Then, to satisfy the antipodal nature of 

the scheme, 8l-2 = it - 2L-1 -• 

• For each m = L — 3, L — 5, ... 2, choose 9m and 8 m -1 = 
n — 9 m from the set 0, given in l |l I [ i. that minimises the sum 
of the condition numbers of the matrices P m and P m_1 , 


• Choose 80 = 0 or #0 = 7r. 











Fig. 2: Plots of the maximum error E max and the mean error _E m ean, 
given in 0 and ( | 1 3| > respectively, for band-limits L < 25. 

Such placement of samples along co-latitude ensures that P™ is 
well-conditioned, resulting in the accurate computation of the pro¬ 
posed SHT as we demonstrate in the next section. 

We have considered the selection of sample positions along co¬ 
latitude from the set of equiangular samples as it has been shown in 
| 20 | | that the selection from the equiangular set results in more ac¬ 
curate computation of the SHT. Furthermore, the equiangular place¬ 
ment of samples along longitude, given in & permits the use of the 
fast Fourier transform for the computation of G m (# n ) 0. We note 
that the proposed sampling scheme can be easily customised for the 
case when the samples along co-latitudes and/or longitude are not 
equiangular. 

Remark 3. The computational complexity of the SHT for the pro¬ 
posed sampling scheme is fO(L 4 )) 120], which is much smaller than 
the complexity of least squares methods (0(L 6 )) used by many ex¬ 
isting sampling schemes mm- 

5. ANALYSIS OF PROPOSED SAMPLING SCHEME 
5.1. Numerical Accuracy 

For the accurate reconstruction of the diffusion signal, using its ex¬ 
pansion in the spherical harmonic basis in 0 - the accuracy of the 
SHT given in 0 is of the upmost importance. A sampling scheme 
is numerically accurate if the inverse SHT of a band-limited signal 
followed by the SHT gives an error between the original and recon¬ 
structed signal on the order of the numerical precision ]18||20| . 

We carry out the following experiment to test the numerical 
accuracy of the proposed sampling scheme and the corresponding 
SHT. We first obtain a band-limited antipodally symmetric test sig¬ 
nal ft in the spectral domain by generating spherical harmonic coef¬ 
ficients (fffT for 0 < £ < L, l even, \m\ < l with real and imag¬ 
inary parts uniformly distributed in the interval [— 1 , 1 ] (( ft)T = 0 
for 0 < £ < L, f!odd, |m| < £). The inverse SHT is then used to 
obtain ft in the spatial domain over the proposed sampling grid (de¬ 
scribed in Section0, followed by the forward SHT to compute the 
spherical harmonic coefficients of the reconstructed signal, (f r )™. 
The experiment is repeated 10 times and the average values for the 
maximum error E mdx and the mean error E mcan , given by 

E max 4 max |(/ t )r-(/ r )ri, (i 2 ) 

E i(/t)r-(/r)n, us) 

£=0 m=—l 



Fig. 3: Plot of the mean error E mean , given in 0 - for the original 
test signal ft and 5 rotated versions of the test signal, and for band- 
limits L < 25. The legend displays the rotation angles in radians in 
the form of Euler angles (a, /3, 7 ), where the rotation is applied to 
the test signal using the zyz rotation convention 0 - Note that the 
reconstruction accuracy is the same for all rotations. 

are recorded and plotted in Fig0for band-limits L < 25 (the band- 
limits of interest in dMRI). It is evident that, although the errors 
_E max and E mcan grow as the band-limit L increases, the errors are on 
the order of numerical precision (less than 10 ~ 14 ), showing that the 
proposed sampling scheme and SHT allows for the accurate recon¬ 
struction of any band-limited antipodally symmetric signal on the 
sphere for band-limits L < 25. 

5.2. Rotational Invariance 

Many existing sampling schemes on the sphere focus on uniform 
sampling of the diffusion signal on the sphere to ensure that the re¬ 
construction accuracy is rotational invariant i.e. does not change sig¬ 
nificantly when the signal or sampling grid is rotated HD- The pro¬ 
posed sampling scheme is not uniform by design, however it does not 
have dense sampling on any region of the sphere. The rotational in¬ 
variance of the reconstruction accuracy of the SHT was tested using 
the numerical accuracy experiment (described in Section [5T| per¬ 
formed on 5 rotated versions of the antipodally symmetric test sig¬ 
nal ft. For each rotation, we randomly obtain Euler angles (a, j3, 7 ) 
from the uniform distributions, where a, 7 G [0, 27r) and fi G [0, tt], 
and apply the rotation using the zyz convention 0 - F 's 0 shows 
the mean reconstruction error E mean , averaged over 10 experiments, 
for the original test signal and 5 rotated versions of the test signal, 
for band-limits L < 25. It can be seen that the reconstruction errors 
are on the same order of magnitude regardless of the rotation angle, 
demonstrating that the accuracy of the scheme does not depend on 
the angle of rotation. 

6. CONCLUSIONS 

In this work we have proposed a new sampling scheme on the sphere 
for the accurate reconstruction of the diffusion signal in dMRI. The 
proposed scheme exploits the antipodal symmetry of the diffusion 
signal and attains an optimal number of samples No = (L + 1)/ 2 , 
given by the degrees of freedom required to represent the diffusion 
signal of band-limit L , in contrast to the existing schemes which re¬ 
quire at least L 2 samples. The smaller number of samples required 
by this scheme will enable a reduction in the scan time required to 
obtain sufficient measurements for the accurate reconstruction of the 
diffusion signal. The proposed sampling scheme enables accurate 


























and rotational invariant computation of the SHT of the diffusion sig¬ 
nal to near machine precision accuracy. Applying the scheme to 
dMRI data and extending the proposed scheme to multiple q-shell 
sampling is the subject of our future work. 
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